function    A204()
format long;

% 统计估计：统计估计（statistical estimation）是统计推断的主要内容，
% 当变量的分布形态未知，根据样本数据对变量的分布形态进行推测（估计）；
% 当变量的分布形式已知，根据样本数据对未知的参数或未知参数的函数进行估计。

% 点估计
% 钜估计
Value = randn(50, 1);
mu_Value = mean(Value);
sigma2_Value = moment(Value, 2);
sigma3_Value = moment(Value, 3);
sigma4_Value = moment(Value, 4);
% 钜估计可视化
subplot(4, 2, 1);
co1_X = 1 : 50;
plot(co1_X, Value, 'b*');
set(gca, 'Fontsize', 12);
str1= {['均值（', num2str(mu_Value), '）'], ...
    ['方差（', num2str(sigma2_Value), '）'], ...
    ['三阶矩（', num2str(sigma3_Value), '）'], ...
    ['四阶矩（', num2str(sigma4_Value), '）']};
text(co1_X(1)+co1_X(50)*0.6, min(Value)+(max(Value)-min(Value))*0.6, str1, 'Fontsize', 12);
title('矩估计');

% 极大似然估计
alpha=0.05;      % alpha为显著性水平，1-alpha为置信水平
x2 = [10 7 13 14 14 9 41 34 37 22 44 33 51 50 37 147 87 23 47 134 83 133 ...
    144 119 103 173 110 130 88 53];
N2 = 30; 
% mu_sigma最大似然估计，mu_sigma_ci为95%置信区间
[mu_sigma, mu_sigma_ci] = mle(x2, 'distribution', 'norm', 'alpha', 0.05); 
% 极大似然估计可视化
subplot(4, 2, 2);
co2_X = 1 : N2;
plot(co2_X, x2, 'b*');
set(gca, 'Fontsize', 12);
str2= {['最大似然估计（', num2str(mu_sigma(1)), '，', num2str(mu_sigma(2)), '）'], ...
    ['置信区间（', num2str(mu_sigma_ci(1, 1)), '，', num2str(mu_sigma_ci(2, 1)), '；'], ...
    ['                ', num2str(mu_sigma_ci(1, 2)), '，', num2str(mu_sigma_ci(2, 2)), '）']};
text(co2_X(1)+co2_X(30)*0.3, min(x2)+(max(x2)-min(x2))*0.6, str2, 'Fontsize', 12);
title('极大似然估计');

% 最小二乘法
x3 =[30 34 39 91 106 73 255 128 248 185 157 135 83 109 95 161 255 260 ...
    302 386 415 550 566 581 857 1032 955 898 708 441];
N3 = 30; 
co3_X = 1 : N3;
[p, S] = polyfit(co3_X, x3, 1);
x3_Fit = polyval(p, co3_X);
% 最小二乘法的可视化
subplot(4, 2, 3);
plot(co3_X, x3, 'b*', co3_X, x3_Fit, 'r-');
set(gca, 'Fontsize', 12);
legend('测量值', ['最小二乘法拟合 灵敏度(', num2str(p(1)), ')，零位输出(', num2str(p(2)), ')']);
title('最小二乘法估计');

% Bayesian估计
total_mu = 3.0;                          % 总体分布概率的均值
total_sigma = 3.6;                       % 总体分布概率的标准差
% 构造100个如下所示的样本：
sample_mu = 2.8;                         % 未知参数分布的样本均值为2.8
sample_sigma = 0.5;                      % 未知参数分布的样本标准差为0.5
sample_num = 100;                        % 样本数为100
sample_X = normrnd(total_mu, total_sigma, 1, sample_num);   % 产生100个样本
sample_muX = mean(sample_X);             % 计算样本均值
Bayes_value = sample_num * sample_sigma^2 * sample_muX/(sample_num * sample_sigma^2 ...
    + total_sigma^2) +total_sigma^2 * sample_mu/(sample_num * sample_sigma^2 ...
    + total_sigma^2);   % Bayesian估计值
rel_error = abs(Bayes_value-total_mu)/total_mu*100;     % 计算Bayesian估计值与样本均值的相对误差
% 计算随着样本逐渐新增，其Bayesian估计值与相对误差值的变化规律
for sample_num = 1: 5000
    new_sample_value = normrnd(total_mu, total_sigma, 1, 1);  % 产生与未知参数分布一样的样本数据
    sample_X = [sample_X, new_sample_value];                  % sample_X逐渐增加样本数据
    sample_muX = mean(sample_X);                              % 样本均值
    Bayes_value(sample_num)=(sample_num+100)*sample_sigma^2*sample_muX ...
        /((sample_num+100)*sample_sigma^2+total_sigma^2+total_sigma^2*sample_mu ...
        /((sample_num+100)*sample_sigma^2+total_sigma^2));   % Bayesian估计值
    % 计算Bayesian估计值与样本均值的相对误差
    rel_error(sample_num) = abs(Bayes_value(sample_num)-total_mu)/total_mu*100; 
end
% 可视化Bayesian估计值与相对误差值随新增样本的变化规律
subplot(4, 2, 4)
sample_num = 1: 5000;
AX = plotyy(sample_num, Bayes_value, sample_num, rel_error);
set(gca, 'Fontsize', 12);
xlabel('增大样本数'); 
set(get(AX(1),'ylabel'),'string', '贝叶斯估计值');
set(get(AX(2),'ylabel'),'string', '相对误差', 'Fontsize', 12);
title('Bayesian估计');

% 区间估计
% 计算随着正态随机分布样本逐渐新增，区间估计的变化规律
muCI = [0];
sigmaCI = [0];
for rnd_num = 1: 90
    rnd_X = normrnd(10, 30, 1, 10+rnd_num);
    [muHat(rnd_num), sigmaHat(rnd_num), per_muCI, per_sigmaCI] = normfit(rnd_X, 0.01);
    muCI_neg(rnd_num) = per_muCI(1);
    muCI_pos(rnd_num) = per_muCI(2);
    sigmaCI_neg(rnd_num) = per_sigmaCI(1);
    sigmaCI_pos(rnd_num) = per_sigmaCI(2);
end
% 可视化区间估计随新增正态随机样本的变化规律
subplot(4, 2, 5)
rnd_num = 1: 90;
errorbar(10+rnd_num, muHat, muCI_neg, muCI_pos, 'bo');
hold on;
errorbar(10+rnd_num, sigmaHat, sigmaCI_neg, sigmaCI_pos, 'r*');
set(gca, 'Fontsize', 12);
xlim([10, 100]);
xlabel('样本数'); 
ylabel({'均值及其区间估计(蓝)', '标准差及其区间估计(红)'});
set(get(AX(2),'ylabel'),'string', '相对误差', 'Fontsize', 12);
title('区间估计');

% 单个正态总体参数的区间估计
% 利用正态逆累积分布函数norminv()
% x = norminv(p,mu,sigma) 基于p中的概率值计算并返回具有均值mu和标准差sigma的正态cdf的逆函数值。
% 计算正态分布参数的最大似然估计MLE，计算对应逆cdf值的置信区间。
rng('default');                   % For reproducibility
n = 1000;                         % Number of samples
x6 = normrnd(5, 2, [n, 1]);       % 从均值为5、标准差为2的正态分布中生成1000个正态随机数。
phat = mle(x6);                   % 使用mle计算分布参数（均值和标准差）的MLE。
muHat = phat(1);
sigmaHat = phat(2);
% 使用normlike估计分布参数的协方差，返回渐近协方差矩阵的逼近。
[~, pCov] = normlike([muHat, sigmaHat], x6);
% 计算在0.5处的逆cdf值及其99%置信区间。like_x6是使用参数为muHat和sigmaHat的正态分布的逆cdf值。
% 使用pCov的muHat和sigmaHat的不确定性，区间[xLo, xUp]是在0.5处计算的逆cdf值的99%置信区间。
[like_x6, x6Lo, x6Up] = norminv(0.5, muHat, sigmaHat, pCov, 0.01);
% 可视化区间估计随新增正态随机样本的变化规律
subplot(4, 2, 6)
num = 0: n-1;
plot(num, x6, 'bo');
hold on;
errorbar(n/2, like_x6, x6Lo, x6Up, 'r*');
legend('测量值', ['逆cdf值(', num2str(like_x6), ')，区间[', num2str(x6Lo), '，', num2str(x6Up), ']']);
set(gca, 'Fontsize', 12);
title('单个正态总体参数的区间估计');

% 核密度估计
% 读取文件A202_data.xlsx中的第一个工作表中的总成绩数据，即F2~F52
score=xlsread('A204_data.xlsx','F2:F52');
score=score(score>0);
%设置核函数分别为Gaussian，Uniform，Triangle和Epanechnikov
%调用ksdensity函数进行四种核密度（正态、均匀、三角和径向基函数）检验
[f_ks1,xi1]=ksdensity(score,'kernel','normal');
[f_ks2,xi2]=ksdensity(score,'kernel','box');
[f_ks3,xi3]=ksdensity(score,'kernel','triangle');
[f_ks4,xi4]=ksdensity(score,'kernel','epanechnikov');
%分布绘制不同核函数所对应的核密度估计图
subplot(4, 2, 7);
plot(xi1,f_ks1,'k','linewidth',2);
hold on;
plot(xi2,f_ks2,'r:','linewidth',2);
plot(xi3,f_ks3,'b-.','linewidth',2);
plot(xi4,f_ks4,'c--','linewidth',2);
xlabel('考试成绩');
ylabel('核密度估计');
legend('正态分布','均匀分布','三角分布','径向函数');
set(gca, 'Fontsize', 12);
title('核密度估计');

% 稳健统计
% 使用Gaussian关联从二元分布生成随机数据点。
rng default
rho = [1, 0.05; 0.05, 1];
u = copularnd('Gaussian', rho, 50);
% 将5个随机选择的观测值修改为异常值。
noise = randperm(50, 5);
u(noise, 1) = u(noise, 1)*10;
% 使用三种方法计算稳健协方差矩阵，即快速最小协方差行列式(MCD), OGK和Olive-Hawkins。
[Sfmcd, Mfmcd, dfmcd, Outfmcd] = robustcov(u);
[Sogk, Mogk, dogk, Outogk] = robustcov(u, 'Method', 'ogk');
[Soh, Moh, doh, Outoh] = robustcov(u, 'Method', 'olivehawkins');
% 利Mahalanobis距离为样本数据计算经典距离值
d_classical = pdist2(u, mean(u), 'mahal');
p = size(u, 2);
chi2quantile = sqrt(chi2inv(0.975, p));
% 为每种稳健协方差计算进行可视化。
subplot(4, 2, 8)
plot(d_classical, dfmcd, 'o')
line([chi2quantile, chi2quantile], [0, 30], 'color', 'r')
line([0, 6], [chi2quantile, chi2quantile], 'color', 'm')
hold on
plot(d_classical(Outfmcd), dfmcd(Outfmcd), 'r+', d_classical(Outogk), dogk(Outogk), 'g*', d_classical(Outoh), doh(Outoh), 'k.')
xlabel('Mahalanobis Distance');
ylabel('Robust Distance');
legend('样本', 'x基线', 'y基线', 'FMCD', 'OGK', 'Olive-Hawkins');
title('计算稳健协方差矩阵的三种方法FMCD、OGK和Olive-Hawkins');
hold off


